Microscopic theory of singlet exciton fission. I. General formulation 

Timothy C. Berkelbach, 1 '0Mark S. Hybertsen, 2 0and David R. Reichman 1 EI 

{) Department of Chemistry, Columbia University, 3000 Broadway, New York, New York 10027, 
USA 

21 Center for Functional Nanomaterials, Brookhaven National Laboratory, Upton, New York 11973-5000, 
USA 

Singlet fission, a spin-allowed energy transfer process generating two triplet excitons from one singlet exciton, has 
the potential to dramatically increase the efficiency of organic solar cells. However, the dynamical mechanism of this 
phenomenon is not fully understood and a complete, microscopic theory of singlet fission is lacking. In this work, we 
assemble the components of a comprehensive microscopic theory of singlet fission that connects excited state quantum 
chemistry calculations with finite-temperature quantum relaxation theory. We elaborate on the distinction between 
localized diabatic and delocalized adiabatic bases for the interpretation of singlet fission experiments in both the time 
and frequency domains. We discuss various approximations to the exact density matrix dynamics and propose Redfield 
theory as an ideal compromise between speed and accuracy for the detailed investigation of singlet fission in dimers, 
clusters, and crystals. Investigations of small model systems based on parameters typical of singlet fission demonstrate 
the numerical accuracy and practical utility of this approach. 



I. INTRODUCTION 

The Shockley-Queisser limit places the maximal efficiency 
of a single-junction solar cell at about 31%.' Promising tech- 
nologies aimed at exceeding this limit include tandem solar 
cells' 3 hot carrier collection, 4 " 6 and multiple exciton gener- 
ation (MEG). 7 ' 8 In MEG and its molecular analogue, singlet 
fission, a single absorbed photon generates two or more ex- 
citons each of lower energy, eventually yielding two or more 
electron-hole pairs. This mechanism results in theoretical so- 
lar cell efficiencies of almost 50% with singlet fission 9 or more 
with MEG. 7 Singlet fission is a particularly promising tech- 
nology in inexpensive organic solar cells, whose efficiencies 
to date remain well below that of their more expensive inor- 
ganic counterparts. Proposals to utilize singlet fission in this 
manner have targeted covalently linked dimers for use in dye- 
sensitized cellsiS as well as crystalline materials for more tra- 
ditional heterojunction cells 

Despite initial reports of singlet fission over 40 years 
agOr^r— an explosion of experimental studies have emerged 
only recently due to the aforementioned potential for pho- 
tovoltaic utility. Singlet fission, as typically measured by 
the observation of triplets in the form of delayed fluores- 
cence (DF) or transient absorption (TA), has been found to 
vary from system to system both in total yield and overall 
timescale, making the search for unifying principles very dif- 
ficult. The authoritative review by Smith and Michl 9 effec- 
tively summarizes the state of the field up to 2010. Since 
then, singlet fission has been further investigated by TA in thin 
films of diphenylisobenzofuran, 16 by DF and TA in crystalline 
tetracener^il^ by time-resolved two-photon photoemission 19 
and TA 20 in crystalline pentacene, by TA in solution and crys- 
talline rubrene, 21 and even by DF and TA in amorphous films 
of diphenyl tetracene. 22 Although still far away from commer- 
cial use in solar cells, singlet fission has been investigated in 
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pentacene-perylene blend films and even successfully incor- 
porated into heterojunction solar cells utilizing phthalocya- 
nine, tetracene, and C6cr^. 

These many enlightening experiments notwithstanding, the 
dynamical mechanism of singlet fission is still not well un- 
derstood. Previous theoretical work has focused almost en- 
tirely on identifying the quantum mechanical states involved, 
including their wavefunction character and energetic order- 
ing. Perhaps most significantly, high-level quantum chemistry 
calculations have identified a multi-exciton state that is com- 
posed of two triplets coupled into an overall spin singleti^i 2 ^ 
The transition to this multi-exciton state is thus spin-allowed 
and should proceed rapidly, while its triplet-triplet character 
suggests that it should naturally relax to separated triplets on 
a longer timescale. Unfortunately, the multi-exciton nature 
of this state implicitly prevents its direct photoexcitation such 
that the state is spectroscopically "dark" and difficult to ob- 
serve. However, time-resolved two-photon photoemission 19 
and transient absorption 18 spectroscopic measurements have 
provided direct and indirect evidence, respectively, of this 
multi-exciton state in ultrafast singlet fission. 

A simple four-electron four-orbital model suggests at least 
two viable mechanisms for the transition from an initially 
excited intramolecular singlet state, Si, to the multi-exciton 
triplet-triplet state, TT^ The first is a mediated mechanism, 
whereby a charge transfer state acts as an intermediate in the 
transition from S \ to TT; theoretical studies of this mecha- 
nism in coupled molecular dimers have considered the static 
electronic parameters^ as well as the real-time dynamics in 
the limit of fast coherent transfer 27 and in the presence of a 
low-frequency solvent bath. 28 Alternatively, a direct mecha- 
nism has also been implicated, whereby the Coulomb poten- 
tial yields a direct interaction between S i and TT, avoiding 
any intermediates. Some authors have invoked such a pro- 
posal to explain fission in crystalline tetracene and pentacene, 
based both on experiment 19 and quantum chemistry calcula- 
tions of clusters^ 

An internally consistent theory of singlet fission phenom- 
ena must comprise a currently nonexistent unification of static 
electronic structure and dynamic relaxation mechanisms. In 
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this first article of a series we pursue this goal, presenting a 
fully microscopic theoretical formalism tailored to the inves- 
tigation of singlet fission in molecular systems. Our goal is the 
identification, extension, and marriage of existing techniques 
of electronic structure theory and microscopic quantum dy- 
namics for the efficient and accurate treatment of singlet fis- 
sion in organic molecules and bulk materials. Such a synthesis 
elucidates experimental results in both the time and frequency 
domains and allows for studies of competing mechanisms as 
well as quantitative predictions. Our approach is related in 
spirit to treatments of excitation energy transfer in photosyn- 
thetic pigment protein complexes, where reduced density ma- 
trix simulations similar to those proposed here have enjoyed 
great success in understanding quantum effects in complex, 
multi-state biological systems. 29-31 The second and third arti- 
cles of this series will make clear the utility of this formalism 
as we investigate singlet fission in dimers and crystals, respec- 
tively. 

The layout of this paper is as follows. In Sec. [II] we present 
a minimal electronic structure model capable of describing all 
relevant states and couplings for the problem of singlet fis- 
sion. We then proceed in Sec.[III]to describe a non-Markovian 
quantum master equation approach for the description of re- 
laxation mechanisms arising from the coupling of electronic 
degrees of freedom with nuclear vibrations in systems under- 
going singlet fission. We present numerical examples of our 
approach in Sec. [IV] as applied to simple model systems of 
singlet fission, benchmarking our results against numerically 
exact calculations. In Sec. [V] we reflect on our approach and 
conclude. 



II. ELECTRONIC STRUCTURE OF SINGLET FISSION 

Our theoretical framework begins with the electronic struc- 
ture of singlet fission chromophore systems in the limit of 
frozen nuclei. In Sec. Ill Al we will emphasize the utility of 
a generically defined diabatic basis of excited states, and in 
particular how they should be interpreted, thereby providing 
a rigorous and important language with which to speak about 
the electronic structure of singlet fission. As a practical, phys- 
ically intuitive example of such a framework, in Sec. Ill Bl we 
will present a limited configuration interaction (CI) descrip- 
tion of these diabatic excited states. This model quantum 
chemical formalism will be recognized as a modest general- 
ization of the picture proposed by Smith and Michl 9 to under- 
stand singlet fission in dimers. Our contributions are to dis- 
tinguish between the diabatic basis and the adiabatic basis, as 
related to both theoretical development and experimental in- 
terpretation, to formalize this procedure within the context of 
ab initio quantum chemistry, and to generalize to the case of 
more than two molecules. We will briefly discuss alternative 
quantum chemistry approaches in Sec. Ill CI before assessing 
the accuracy and summarizing the proposed electronic struc- 
ture formalism in Sec. Ill Dl 



A. The basis of diabatic states 

Consider M molecules with a total of N electrons. We start 
by defining a basis of diabatic electronic states, i.e. those 
states whose quantum mechanical character is well-defined 
and presumed to be independent of the molecular geometry. 
For an accessible discussion of diabatic states in the context 
of electron transfer, see Ref . 

The desired basis begins with the exact many-electron 
wavefunction corresponding to the (singlet) ground state of 
the system, *F,y (ri,r2, . . . , r^) = (r|5o). All remaining states 
in the minimal singlet fission diabatic basis are excitations 
above the ground state, although they are not in general eigen- 
states of the electronic Hamiltonian. Specifically, the basis 
must include all M excited states which may be character- 
ized by molecule m being in its first excited singlet (or S i ) 
state, \(S i),„), sometimes referred to as Frenkel excitations. 
The next class of excited states are charge-transfer (CT) states, 
where molecule m has a single positive charge and molecule n 
has a single negative charge, denoted |C„,A„) (C for cation and 
A for anion). The final class of necessary excited states are 
those multi-exciton states described as being a spin-adapted 
combination of triplet excitations on molecules m and n, form- 
ing an overall singlet, \T m T n ). 

Double excitations that instead couple two singlets are pre- 
sumed too high in energy to be relevant for singlet fission and 
are consequently neglected, as are all double excitations in- 
volving more than two molecules. States of differing multi- 
plicity (such as triplet- and quintet-coupled triplets, 3 7T and 
5 TT) can also be included. While a purely electronic Hamil- 
tonian does not couple such states of different multiplicity, the 
spin dipole-dipole Hamiltonian does. 9 We neglect this Hamil- 
tonian in our present formalism because it is significantly 
weaker, and we instead only focus on the short-time forma- 
tion of the spin-singlet multi-exciton state, \T m T„). However, 
including spin dipole-dipole and Zeeman interactions in the 
presence of a magnetic field would be important but straight- 
forward modifications to study such long-time effects. 

Although these many-electron basis states, which we will 
denote generically by |/) and are not eigenstates of the 
electronic Hamiltonian and do not constitute a complete basis, 
one may consider the projection of the true electronic Hamil- 
tonian onto this basis, 

H el *J]\i}(i\H el \j)(j\- (1) 

ij 

We wish to emphasize that the states defined above merely 
constitute a physically-motivated (diabatic) basis. Many ex- 
perimental measurements instead observe the adiabatic basis, 
which is is approximately obtained from the diagonalization 
of the above electronic Hamiltonian, H e /, yielding eigenstates 
that are a mixture of the diabatic states. 33 Therefore, one must 
exercise great caution when speaking about the character of 
observed states, e.g. "charge-transfer," and when discussing 
results and proposing mechanisms in terms of these states. 
This distinction must also be kept in mind for traditional ex- 
cited state electronic structure calculations, which inherently 
probe the adiabatic, and not the diabatic basis. 
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In light of the above proviso, one may naturally question 
the utility of this diabatic basis. We propose three reasons to 
begin the theoretical development of singlet fission from this 
basis: 

i. ) The physical character of the basis states aids in in- 
terpreting the nature of observed eigenstates, allowing for a 
means to quantify statements such as "a mixture of charge- 
transfer and Frenkel excitations." This latter example will 
play a prominent role in our future work on singlet fission 
in crystals. Similarly, this principle underlies the coherent su- 
perposition approximation recently proposed to explain MEG 
in nanocrystals 34,35 and singlet fission in pentacene, 19 wherein 
single- and multi-exciton (diabatic) states are coupled to yield 
an eigenstate that is a superposition of the two. 

ii. ) The local diabatic basis can yield accurate results which 
computationally scale very favorably. Proximity arguments 
alone can naturally suggest coupling terms that may be ap- 
proximated or neglected entirely. Using such approximations, 
to be discussed in more detail in our next paper, one may eas- 
ily build up a large molecular aggregate Hamiltonian using 
only diabatic energies and couplings from monomers, dimers, 
or small clusters, which may be computed with very high ac- 
curacy. This philosophy is reminiscent of fragmentation meth- 
ods in the pursuit of linear scaling quantum chemistry. 36 

iii. ) Lastly, the molecular character of the diabatic ba- 
sis allows for a straightforward extension to include cou- 
pling to molecular vibrations, which naturally separate into 
intramolecular and intermolecular modes, as we detail in the 
following section. 

Clearly, the accurate construction of diabatic states marks 
an important research goal for ab initio simulations of sin- 
glet fission. While our approach here and henceforth em- 
ploys a constructive strategy, i.e. a direct construction of 
diabats without explicit reference to the adiabatic states of 
the extended system, an alternative route would employ de- 
ductive strategies that attempt to obtain approximate diabats 
given a set of adiabats. This latter set of states is more easily 
obtained at high accuracy from existing quantum chemistry 
methods, although the non-uniqueness of this diabatization 
procedure results in various competing methods with subtle 
differences . 32,37 ' 38 In any case, the framework presented here 
is not limited to the Cl-type model Hamiltonian outlined be- 
low, and more accurate diabatic states, as might be obtained 
from multi-reference quantum chemistry methods, can be nat- 
urally incorporated into the dynamical scheme to be discussed 
in Sec. HID 

B. A minimal, truncated CI basis 

The accurate quantum mechanical calculation of excited 
states in large molecular systems is still a difficult challenge 
(see Refs.O and HI for examples of recent high-level quan- 
tum chemistry calculations as applied to singlet fission) and 
thus we consider here the simplest possible model Hamilto- 
nian approach that captures the essential physics contained 
in the diabatic framework outlined above. Specifically, we 



consider the minimal active space of all Hartree-Fock (HF), 
or HF-like, highest occupied and lowest unoccupied molecu- 
lar orbitals (HOMOs and LUMOs) of the isolated molecules; 
extension to include additional frontier orbitals is straightfor- 
ward. We furthermore restrict the electronic structure calcula- 
tion to all single and select double excitations, the latter ensur- 
ing treatment of the bi-excitonic triplet-triplet state. If done as 
a purely ab initio theory, this approach would be somewhat 
akin to configuration interaction 39 with single and (select) 
double excitations (CISD) with the frozen core and deleted 
virtuals approximations, or alternatively a type of (severely) 
restricted active space CISD. However, our formalism differs 
slightly in that we consider excitations among the isolated 
molecular orbitals, rather than among the HF orbitals of the 
full interacting system. 

To make our description more precise, we define the cre- 
ation (annihilation) operator for the HOMO of molecule m 
with spin cr as c* H {cH,m,a) and likewise for the LUMO. The 
ground state is thus taken to be 

M 

iso> = n n o°> 

m=\ tr=T,l 

where |0) is the vacuum state of inactive core orbitals, thus fill- 
ing the HOMO of all molecules. As discussed above, this state 
is not the result of a self-consistent HF procedure. From this 
ground state, we will generate the three types of excited states 
described above in Sec. Ill Al Because the electronic Hamil- 
tonian is spin-conserving, we take symmetry-adapted linear 
combinations of select excitations to generate simultaneous 
eigenstates of the S z and S 2 operators, with eigenvalues of 
for both (sometimes called configuration state functions). 

The first type of state is the local Frenkel singlet excitation 
on molecule m, given by 

\(S i) m > = — (c[„ h fH,m,t + c i, m ,X C «."'.i) I 5 0>- (3) 

In addition to the above intramolecular excitation, the single 
excitations also generate our second type of state, namely the 
intermolecular charge-transfer excitation obtained by exciting 
an electron from the HOMO of molecule m to the LUMO of 
molecule n{n + m), 

\C m A n ) = — (c[ n fH,m,T + C [^i C H,m,l) \S 0>, (4) 

where C and A denote the cationic and anionic species, respec- 
tively. The above two types of excited states combine to yield 
all possible single excitations, so that stopping at this point 
would constitute a full Cl-singles (CIS) within the HOMO- 
LUMO space. 

However as discussed above, the problem of singlet fission 
necessarily requires our third type of state, a double excita- 
tion coupling two intramolecular triplet excitations into a state 
with overall singlet character, 
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TmTn) = ^[ 2C i.»4 C L,T Cff -"-t C ^i +2c L,t C lm,i Cff ." 



|CH, m , T - C L m i C L n i C HM C H ,m,l 



+ C, ,Cr 

L,m,l L. 



(5) 



Finally, we point out that because the molecular orbitals 
of distinct isolated molecules are not necessarily orthogonal 
to one another, the use of creation and annihilation operators 
acting in the space of these orbitals is not strictly rigorous. 
While one could imagine employing suitably orthogonalized 
molecular orbitals that retain the localized nature of isolated 
orbitals, the actual overlap in molecular dimers and crystals is 
often negligibly small, thus justifying the theory in its present 
form. 

Having defined a set of diabatic basis states, it thus remains 
to calculate all matrix elements of the electronic Hamiltonian, 
(i\H e i\j). While the calculation is straightforward, the results 
are cumbersome, and so we include the explicit results in 
App. [A] As discussed more below in Sec. Ill Dl the diagonal 
matrix elements (energies) are only expected to be of qualita- 
tive accuracy but can provide useful insight, and likewise for 
the off-diagonal elements (couplings). For example, the cou- 
plings naturally separate into two classes: those containing 
one-electron integrals and those containing only two-electron 
integrals. The one-electron integrals include the simple ki- 
netic energy term, describing favorable charge derealization, 
or "hopping." Such one-electron integrals are expected to be 
one or more orders of magnitude larger than the two-electron 
ones. Reasonable estimates for typical singlet fission chro- 
mophores in close proximity are 50-100 meV for one-electron 
integrals and 5 meV or less for two-electron integrals. These 
simple analytical expressions and order of magnitude esti- 
mates contribute to the interpretation of singlet fission in terms 
of mediated and direct mechanisms. Qualitatively, the me- 
diated mechanism proceeds via two one-electron processes, 
whereas the direct mechanism proceeds by one two-electron 
process. Which of these two diametric mechanisms prevails 
in a given system of interest will depend sensitively on the 
relative energies of the diabatic states and the dynamics of the 
nuclear degrees of freedom, as will be discussed in Sec.lIVI 



C. Aside regarding wavefunction free methods 

Although the formalism here has employed the HF orbitals 
to construct a many-electron basis, we pause to consider some 
alternatives. At least at the level of single excitations, many 
other electronic structure theories can be reduced to an eigen- 
value equation for the transition energies, much like CIS. Note 
that the CIS theory amounts to the diagonalization of an effec- 
tive two-particle Hamiltonian, 

H m = 5 ik 5j, (sj - £,) + [ft - fj) % jM , (6) 

where /; is the ground-state occupancy of orbital i and 

% jM = 2{il\jk) - {il\W(r^r 2 )\kj). (7) 



Clearly, for CIS, W(ri,r 2 ) = r\\. Physically, the vertex 7C 
describes the interaction between single-particle excitations 
i — > j and k — > I. If the original single-particle states are 
not a good approximation to the quasiparticles of the system, 
as determined e.g. by comparison with electron affinity and 
ionization energies, then the HF excitations are in some sense 
a poor starting point on which to build interactions. In other 
words, the true many-body excitations will require contribu- 
tions from many single-particle excitations. 

Instead, one could start from a ground state density func- 
tional theory (DFT) calculation of the isolated molecules, and 
then consider excitations within the Kohn-Sham (KS) orbitals; 
we will not dwell here on the physical reasons for which the 
KS orbitals may be better single particle states. Suffice it to 
say that this approach is adopted in time-dependent DFT (TD- 
DFT) 40-42 and many-body Green's function approaches'^ 
both of which typically yield results superior to those of 
HF-based CIS, finding many-body excitations strongly dom- 
inated by far fewer single-particle excitations. For example, 
the Green's function based Bethe-Salpeter equation (BSE) in 
practice adopts the form Eq. <JSJ with perturbatively corrected 
orbital energies and a statically screened interaction, 

W( ri ,r 2 )~ ^ dre-\r Xl r 1 cj = Q)\r-r 2 r\ (8) 

where e(ri,r2,w) is the frequency-dependent dielectric 
function. 44 In the crude limit where e _1 (ri,r2,w = 0) = 
e~ l 5(r\ -r 2 ), with e a dielectric constant, one arrives at simply 
static screening of the direct Coulomb term, 

<K ijM = 2(il\jk) - e-\il\kj). (9) 

On the other hand, the inclusion of doubly excited states 
presents an ongoing challenge to these former methodolo- 
gies, making them difficult to employ in an internally con- 
sistent theory of singlet fission. Recent work has demon- 
strated that higher excitations only arise in the above theo- 
ries with the retention of a frequency dependent interaction 
kernel. 46,47 Specifically, one solves the eigenvalue-like equa- 
tion, H((jS)c = (DC, with 

Hij,ki(u>) = S ik 6 ji [ej - e) 

(10) 

+ (fi - fj) [2(il\jk) - (U\W(r u r 2 , u)\kj)\ . 

The operator WicS) is related to the dynamically screened di- 
electric function in BSE and to the exchange-correlation ker- 
nel in TD-DFT (i.e. the adiabatic approximation precludes 
observation of multiple excitations in TD-DFT). We consider 
this a very interesting research focus for singlet fission and a 
subject of future work. 
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D. Accuracy and summary 



A. System-bath Hamiltonian 



Based on the preceding discussion, the numerical results 
of the approach proposed in Sec. MB I will only be of qual- 
itative accuracy. The diagonal energies, (i\H e i\i), should be 
considered estimates of their true values, with some leeway 
for semi-empirical adjustment. For example, we find that the 
gas-phase So — ♦ Si transition energy of pentacene predicted 
by the above approach with a 6-31G(d) basis set is 3.76 eV, to 
be compared to the experimental value of approximately 2.3 
eV i 48 i 49 Including the additional dynamical correlation arising 
from the frozen orbitals (i.e. not just the HOMO and LUMO) 
yields the improved value of 2.81 eV. While TD-DFT is typ- 
ically expected to be an improvement, it was shown previ- 
ously to predict values of 1.64 eV and 1.90 eV, for the PBE 
and B3LYP functionals, respectively. 50 Thus, even purport- 
edly sophisticated methods yield excitation energies with er- 
rors ranging from 0.4 to 0.7 eV>^2' 51 Interestingly, the ad-hoc 
BSE-like prescription, Eq. ©, using only the DFT HOMO 
and LUMO from B3LYP and the dielectric constant of pen- 
tacene e = 3.6, predicts a transition energy of 2.95 eV, much 
improved from the HOMO-LUMO CIS result (note that there 
is, however, no a priori reason that the dielectric constant for 
bulk pentacene should be physically meaningful for a single 
molecule). Only multi-reference perturbation theory^ 4 - and 
full many-body GW/BSE calculations 52 yield quantitative ac- 
curacy, predicting 2.1 and 2.2 eV, respectively. Similarly, the 
electronic couplings in this basis may not be quantitatively ac- 
curate, but have already been shown in other work to provide 
useful qualitative insight into the efficiency of singlet fission 
through investigation of their relative magnitudes 26 and de- 
pendence on molecular orientation. 9 It may thus be permissi- 
ble to uniformly scale the electronic coupling matrix elements 
when investigating singlet fission. 

To summarize, we argue that the diabatic basis, comprising 
states that are easily characterized and energies and couplings 
that are straightforwardly calculated, acts as the crucial con- 
ceptual intermediate between high-level quantum chemistry 
calculations, which inherently yield electronically adiabatic 
states that are difficult to characterize, and microscopic quan- 
tum master equations, which are required to accurately treat 
thermally induced relaxation effects, the topic of the next sec- 
tion. 



SYSTEM-BATH QUANTUM DYNAMICS 



In this section, we consider the coupling of electronic (sys- 
tem) and nuclear (bath) degrees of freedom. Although the 
treatment is relatively standard and can be found in textbooks, 
see e.g. Ref. 53, we include the derivation in App.lBl to em- 
phasize the microscopic connection to the diabatic basis in- 
troduced above. The result is the system-bath Hamiltonian 
described in Sec. MI Al 



To include the effects of electron-phonon coupling, we em- 
ploy a system-bath type Hamiltonian 



H t ot = H e l + Hel-ph + H p i, 



(11) 



with the electronic Hamiltonian calculated at the ground-state 
geometry in terms of the diabatic states described in Sec. Ml 



H el = ^ \i)Ei(i\ + J] lOVyOI, 



(12) 



the bilinear electron-phonon coupling, 



Hel-ph - Ckjflk (13) 

i k ij k 



with Ckfl = 0, and the free phonon Hamiltonian, 

^2 



H. 



ph 



z 



Pk A 1 2.2 

y + 2"** 



(14) 



In the above, i and j index the diabatic electronic basis states, 
and k indexes both the inter- and intra-molecular (ground 
state) normal modes of the system. 

The molecular vibrations and phonons are completely de- 
scribed by their spectral density, 



7T v-" 1 C k ij 
Jij{o)) = - > . —t-oia) - OJk) 



(15) 



Physically, the spectral density encodes the distribution of 
normal mode frequencies weighted by the strength with which 
each mode couples to the energy level of diabatic state i 
(Ju(u))) or to the electronic coupling between states i and j 
{Jjj(u))). In practice, spectral densities (obtained in a manner 
to be described) are usually fit to a numerically convenient 
functional form, J(a>) = AF((l>/Q.), parametrized by the reor- 
ganization energy, A - n~ x J d<jL>J(a))lu and a characteristic 
frequency O. 

Atomistically, the spectral densities may be calculated 
through a combination of classical molecular mechanics and 
quantum chemistry calculations. In one approach, a direct 
diagonalization of the molecular mechanics Hessian yields 
phonon frequencies a>k and displacement vectors, and quan- 
tum chemistry calculations along these displacements pro- 
duce the coupling constants c*. Such an approach has been 
adopted recently by Girlando et al. in studies of electron and 
hole transport in rubrene 54 and pentacene 55 crystals. 

Alternatively, by appealing to the quantum-classical cor- 
respondence of harmonic oscillators, which are presumed to 
compose the nuclear bath, one may show that the spectral den- 
sity can be obtained from the Fourier cosine transform of a 
classical correlation function, 56,57 



Jij{aS) 



o 



Cfj(t) cos(wf) 



(16) 
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where C?/(f) = (SEi(t)6Ei(Q))j is the energy gap fluctu- 
ation correlation function and C c £ (f) = (6Vij(t)6Vij(0))j 
is the electronic coupling fluctuation correlation function 
(6X = X - (X)j). This latter approach has been exten- 
sively pursued in the present context of organic materials by 
Troisi and co-workers, who have focused on the fluctuations 
and spectral properties of the electronic coupling in DNA, 58 
pentacene crystals, 59 and the discotic phase of hexabenzo- 
coronene derivatives.— 

Given the expense of accurate ab initio quantum chemistry 
methods, frequent calculations along the course of a molec- 
ular dynamics trajectory are clearly prohibitive. As such, 
it is common to adopt a semi-empirical quantum chemical 
method, such as the spectroscopic parametrization of INDO 
(intermediate neglect of differential overlap), which has an 
impressive accuracy to cost ratio allowing for the rapid col- 
lection of sufficient statistics. While one could in principle 
calculate all the diabatic matrix elements defined in App. [A] 
we note that the diagonal elements are dominated by the bare 
orbital energies and the off-diagonal coupling matrix elements 
are dominated by the one-electron coupling. Thus it is reason- 
able to assume that the stochastic properties (fluctuation mag- 
nitude and timescale) of the full matrix element are equivalent 
to those of its one-electron terms. These latter properties are 
more commonly evaluated in the literature, due to their role in 
electron and hole transport of organic materials. 

The last topic of discussion concerns the correlation of dif- 
ferent bath modes, for example the extent to which the fluc- 
tuations of the diabatic energy of state i are correlated with 
those of state j. Although there is surely some degree of cor- 
relation, positive or negative, the effect of its inclusion on the 
subsequent dynamics is debatable. In particular, while some 
studies have attempted to implicate correlated bath modes in 
efficient biological energy transport, 61,62 molecular dynamics 
simulations of photosynthetic complexes show no significant 
correlations. 63,64 Lacking any firm evidence either way for the 
problem of singlet fission, we will let the correlation of differ- 
ent bath modes be dictated by numerical convenience (usually 
preferring the completely uncorrelated scenario), though it is 
a topic worthy of further investigation. 



B. Reduced density matrix dynamics 

The dynamics of the coupled electron system and phonon 
bath is given by the Liouville-von Neumann equation for the 
total density matrix, W(t), 



dW(t) 
dt 



= -i [H lot , W{t)\ , 



(17) 



the exact solution of which is prohibitively difficult due to 
the large Hilbert space associated with the phonon degrees of 
freedom. However, as long as one is only interested in elec- 
tronic observables, great simplification occurs when consid- 
ering the reduced density matrix (RDM) of the system, p(f), 
obtained by averaging the total density matrix over the phonon 
degrees of freedom, i.e. p(f) - Tr p /,{VK(f)}. The diagonal ele- 
ments of this matrix, p,,(f) = (i\p(t)\i), are the populations of 



state i and the off-diagonal elements, p,y(f) = (i\p(f)\j), are the 
coherences between states i and j. 

A variety of methods exist for the determination of the 
RDM, each with its own caveats. Although impressive 
progress has been made in the development of numerically 
exact methods - including path-integral techniques^— the 
multi-configurational time-dependent Hartree ansatzr 9 -^ and 
hierarchical equations of motion 29,72 - we will limit ourselves 
here to approximate methods which are more physically trans- 
parent and more readily applied to very large systems, as will 
be demonstrated in our future work on clusters and crystals. 

Approximate methods are generally perturbative in nature, 
and differ in their choice of perturbative parameter. Clearly, 
the physical problem at hand should dictate the appropriate 
small parameter, thus controlling the accuracy of the perturba- 
tive approximation. The first common approach is to treat the 
electronic couplings in the diabatic basis, V, ; , to second order 
in perturbation theory, while treating the system-bath interac- 
tion exactly; this philosophy comprises Marcus 73 and Forster- 
Dexter 74,75 theories, as well as the more sophisticated nonin- 
teracting blip approximation^ 6 * 7 - 7 - Although this methodology 
has been previously employed in a study of CT-mediated sin- 
glet fission, 28 to be discussed later in this paper, we will ad- 
vocate for an alternative approach which treats the electronic 
couplings exactly in exchange for a perturbative treatment of 
the system-bath interaction. The relative merits of the two ap- 
proaches will be contrasted in Sec. lIVI 

Specifically, we shall pursue the use of a Redfield-like 
equation^— in either its non-Markovian or Markovian form. 
Non-Markovian prescriptions can either take a time-local or 
time-nonlocal form, which corresponds to a series resumma- 
tion in terms of different time-ordered cumulants.^— We will 
present equations for the time-local form (or partial ordering 
prescription), though treatment of singlet fission dynamics in 
terms of the alternative time-nonlocal form (or complete or- 
dering prescription) would be straightforward.— 

In the basis of electronic eigenstates, H e i\a) = %u) a \a), and 
adopting the notation u> a p — a> a - o)p, the time-local Redfield 
equation is given by 



dp aj3 (t) 
dt 



— —iteapPapit) + Rgpyd(t)Pyd(t), 



(18) 



where the initial condition of the total density matrix implic- 
itly takes the factorized form W(0) = p(0)e~ 6 " hIkBT /Z ph , with 
the phonon partition function Z p h = Tr p h{e~ Hp ^ kBT }. This ini- 
tial condition is consistent with an impulsive Franck-Condon 
excitation at time t = 0. In Eq. (fT8l . the first term on the 
right hand side is responsible for coherent energy transfer 
whereas the second term is responsible for population relax- 
ation, coherence transfer and dephasing, and more compli- 
cated population-to-coherence transfer processes. Explicit ex- 
pressions for the Redfield tensor elements, which include in- 
tegrals over thermal bath correlation functions, can be found 
in App.0 

In the limit where the bath relaxation takes place signif- 
icantly faster than that of the electronic system, the time- 
dependent Redfield equation is well approximated by its 
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Markovian form, obtained from Eq. (fT8T l by the replacement 
Ra/3ys(t) - > Rafiydi 00 )- This approximation clearly simplifies 
the form of the density matrix equation and provides a direct 
microscopic route to dephasing and relaxation rates which 
are often employed in other contexts as phenomenological pa- 
rameters. This Markovian approximation should be carefully 
checked for its accuracy in each situation of interest. 

As is commonly done in theories of exciton transport, one 
may furthermore employ the secular approximation to the 
Markovian Redfield equation, which preserves the positivity 
of the RDM, i.e. p„ > O i 53 i 79,81 ' 82 The secular approximation 
amounts to neglecting those elements of the Redfield tensor, 
Ra/iys, for which \cj a p - o) y g\ + 0. In doing so, one decou- 
ples the dynamical evolution of populations and coherences 
in the eigenstate basis. In addition to preserving positivity, the 
secular approximation furthermore guarantees that the system 
RDM approaches thermal equilibrium at long times, i.e. 

p{t -> co) = e- 6 °> /k » T /Z e i, (19) 

which is the correct physical result outside regimes of strong 
system-bath coupling. 

IV. APPLICABILITY AND ACCURACY OF REDFIELD 
THEORY FOR SINGLET FISSION DYNAMICS 

Although the presentation of system-bath dynamics up to 
this point has been largely generic, we now thoroughly dis- 
cuss the applicability of the Redfield equation to the specific 
problem of singlet fission in organic systems. We must first 
acknowledge the potential disadvantages of the Redfield treat- 
ment and the extent to which they affect the reliability of such 
calculations. The main approximation inherent in this ap- 
proach is the assumption of weak coupling between the elec- 
tronic and vibrational degrees of freedom. This coupling can 
be quantified approximately by the ratio of the magnitude of 
fluctuations in the nuclei, Ah, to the frequency of these fluctu- 
ations, fiy, If the dimensionless ratio Aij/hCln is small, then 
the Redfield approximation should be a good one. Whether 
this inequality holds or not will depend on the specific system 
under study. 

As a prototypical singlet fission material, consider pen- 
tacene, which we will study in the follow papers. The di- 
agonal reorganization energy and frequencies of the electron- 
phonon coupling have been calculated by quantum chemical 
and molecular dynamics methods 55 to be approximately 50 
meV and 170 meV, respectively; note that this latter value cor- 
responds to the well-known m 1400 cm" 1 aromatic stretching 
mode. Thus we see that the ratio of the two is indeed sig- 
nificantly smaller than one, and the Redfield equation should 
be reasonably accurate. As a general rule, smaller molecules 
will undergo larger geometry distortions in excited states, i.e. 
larger An, and therefore the Redfield approach may break 
down. 

The Markov approximation to the Redfield equation, as dis- 
cussed above, relies on timescale separation between elec- 
tronic and nuclear relaxation, and thus one must compare the 
electronic frequencies ojh to those of the vibrations Qy. With 



electronic frequencies on the order of 50 meV, and again vi- 
brational frequencies of 170 meV, even the Markov approxi- 
mation should be reasonably reliable. In addition to this math- 
ematical argument, there is also a more physical implication 
of the Markov approximation: although the time-dependent 
variants of the Redfield equation can describe multi-phonon 
effects to varying degrees of accuracy, the Markov approx- 
imation inherently describes only single-phonon relaxation 
mechanisms. This deficiency can be readily seen in the adia- 
batic population relaxation rate, R aa /3p, which is proportional 
to Jij((o a p), so that all transition frequencies must be matched 
by a single phonon frequency in the spectral density. 

Potential pitfalls behind us, we now enumerate the many 
advantages of the Redfield formalism. The first obvious ad- 
vantage, which is shared by a variety of other perturbative 
methods, is the clear microscopic formalism. While density 
matrix calculations have been employed for theoretical stud- 
ies of MEG 35 and singlet fission— as well as for fitting exper- 
imental singlet fission data, 19 such dynamical investigations 
have been essentially phenomenological to date. In the ap- 
proach advocated here, the electronic structure methodology 
is directly connected to the molecular structure and micro- 
scopic relaxation mechanisms. The Redfield tensor prescribes 
temperature-dependent population relaxation and coherence 
dephasing rates which can be traced back to the physical vi- 
brations of the system under study. When necessary, the time- 
dependent Redfield variants even yield non-Markovian behav- 
ior which of course cannot be captured with phenomenologi- 
cal time-independent rates. 

In addition, like all master equation methods, the Redfield 
approach scales very favorably in a computational sense. The 
additional adoption of the Markov and secular approximations 
further reduces the computational cost. Needless to say, none 
of the numerically exact methods alluded to above takes on a 
simple master equation form and thus each has a significant 
computational overhead with a scaling that depends on the 
details of the method. 

The theoretical study most similar in spirit to our own is 
that of Teichen and Eaves 28 who sought to quantify the ef- 
fects of a generally non-Markovian bath of low-frequency sol- 
vent degrees of freedom and its implications for CT-mediated 
singlet fission. These authors employed methodology similar 
to the noninteracting blip approximation (NIBA) known from 
spin-boson theory ^£01 previously generalized to the case of 
multilevel system a 66 i 87 and recently extended to situations of 
slow, near-classical bath modes 88,89 in a time-nonlocal formal- 
ism. The time-nonlocal methodology, henceforth referred to 
as NIBA even for multilevel systems, yields a non-Markovian 
master equation for the populations of the RDM in the dia- 
batic basis, P,(f) = pa(t), 

= 2 X dsKij(f ' s)Pj{s) ' (20) 

where 

K u (t, s)=2\ Vj/Re (exp {-iHft) exp (i&fis))^ , (21) 
and// 7 "' = (i\H""\i). Teichen and Eaves instead considered the 
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time-local version of this theory, P,(f) = I t {t) + XjRij(i)Pj(i), 
but the two methods should give similar results, and are iden- 
tical in the Markovian limit. 

As alluded to previously, the NIBA-type master equations 
are perturbative in the electronic couplings, Vy, and thus the 
diabatic basis is in some sense a preferred basis. The nonper- 
turbative effects of strong electronic coupling, yielding signif- 
icant mixing in the adiabatic basis, cannot be described by the 
NIBA theory. Accordingly, as a theory for populations only, 
NIBA makes no prediction about coherence variables, pij(t), 
preventing the transformation to any other electronic basis. As 
described in more detail in App. [D] spectroscopy probes the 
dynamics of coherences in the adiabatic basis, and as such is 
completely beyond reach of NIBA-based theories. 

On the contrary, the nonperturbative nature of Redfield the- 
ory with respect to the electronic Hamiltonian allows for an 
exact solution of the electronic structure problem in exchange 
for an approximate treatment of the system-bath interaction. 
Thus all questions concerning derealization, quantum coher- 
ence, and spectroscopy are readily addressed with the Red- 
field framework, as long as the system-bath coupling is not 
too large. Even in regimes where the time-dependent Red- 
field theory is pushed past its limits of validity, the secu- 
lar and Markovian approximations yield a numerically stable 
Lindblad-type master equation, with microscopically-derived 
relaxation and coherence dephasing rates. Interesting recent 
work has formulated a stable theory which reinserts micro- 
scopic expressions for the population and coherence coupling 
within the Lindblad formalism. 90 



A. Results for population dynamics 

Given the advantages of a Redfield-type approach with re- 
spect to the flexibility of treating populations and coherences 
on equal footing in either the diabatic or adiabatic bases, as 
well as the ability to treat extremely large systems, it is nat- 
ural to ask if such an approach is accurate for typical singlet 
fission systems of current interest. Here we show with small 
model systems that indeed treating the system-bath coupling 
as a perturbative parameter should yield semi-quantitative ac- 
curacy over a wide range of scenarios rooted physically in 
the expected parameter space of acene systems. In all of 
the following results on diabatic population dynamics, we 
make comparison with the numerically exact but computa- 
tionally expensive hierarchical equations of motion (HEOM) 
methodology ^ 29 ' 72,91 as implemented in the Parallel Hierarchy 
Integrator (PHI). 92 To achieve convergence, we truncated the 
hierarchy at L = 5 and required K — 3 terms in the Matsubara 
expansion. 

We begin with a two-state system, which in the context 
of singlet fission may be taken as a model for the direct, 
Coulomb-mediated fission mechanism. The first state is the 
photoexcited initial singlet, Si, and the second state is the 
multi-exciton configuration, TT. The initial condition is 
p(0) = [Si)(Si| and the dynamics proceeds based on the pa- 
rameters of the system-bath Hamiltonian defined above. The 
system-bath coupling will be chosen to take the simple form 



(a) A = 25 meV (b) A = 50 meV 
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(c)A=100meV (d)A = 200meV 




Time t [fs] Time t [Is] 

FIG. 1 . Population dynamics of the two-state singlet fission model described 
in the text with £y, - E TT = 75 meV, V = 50 meV, HQ. = 150 meV (fl -1 ss 
4 fs), and T = 300 K (kgT as 26 meV), for increasing system-bath coupling 
strength. 



Hei-ph = Zi=Si,Tr'Lk,i\i)ck,iqk,i( i l i- e - linear, diagonal cou- 
pling to uncorrected bath degrees of freedom. The baths will 
be characterized by identical Ohmic spectral densities with a 
Lorentzian cutoff (sometimes referred to as the overdamped 
Brownian oscillator model), Ju(co) - 2AQ.a>/(a> 2 + Q 2 ). 

For concreteness, we will fix the fission to be mildly 
exothermic, Es, - Ejj = 75 meV, with a bath cutoff fre- 
quency HQ. = 150 meV (characteristic of aromatic molecules) 
and temperature T = 300 K (IcbT w 26 meV). However, to in- 
vestigate the perturbative accuracy of the Redfield and NIBA 
equations, we will scan the reorganization energy, A, and elec- 
tronic coupling, V. 

In Fig.Q] we plot the singlet population dynamics for a va- 
riety of reorganization energies, with the electronic coupling 
fixed at V = 50 meV. When the system-bath coupling is weak, 
the quantum beating is dominant and overall relaxation is 
slow. It is clear that the timescale of beating should not be con- 
fused with the relaxation timescale; only the former is accessi- 
ble within static electronic structure calculations, whereas the 
latter requires explicit treatment of the vibrational degrees of 
freedom. As the coupling is increased, all theories correctly 
predict that the oscillations become damped and the relaxation 
to TT proceeds more quickly. We see that as the system-bath 
interaction becomes large, the time-local Redfield result be- 
comes inaccurate at long times, even leading to unphysical 
negative populations. However, to some extent, the secular 
and Markov approximations to the Redfield equation prevent 
such a catastrophe, leading to much more reasonable equilib- 
rium populations. The non-Markovian behavior, on the other 
hand, can be observed in the short-time dynamics, which are 
always correctly described by the time-local Redfield equa- 
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(a) V = 1 meV (b) V = 5 meV 
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FIG. 2. The same as in Fig. [T] but with A = 25 meV and scanning the elec- 
tronic coupling, V. Note the changing scale of the time axis. 



tion, but not by its secular, Markovian counterpart. In con- 
trast to the breakdown behavior of the Redfield equation, the 
NIBA dynamics retain their relative accuracy at all values of 
the system-bath coupling. This result is to be expected in as 
much as the NIBA theory treats the system-bath interaction 
exactly. Rather, the NIBA theory is perturbative in the elec- 
tronic coupling, which is unchanged in all panels of Fig. \T\ 
However, NIBA can be seen to consistently underestimate the 
equilibrium population of S \ . This tendency towards extreme 
localization in biased systems is a known deficiency of meth- 
ods that are perturbative in the electronic coupling i 76 ' 77 i 93 

The rapid singlet fission observed in Fig. Q] with 100 fs 
timescales, is due to the rather large value of the direct cou- 
pling matrix element, V. As alluded to previously, this number 
is likely significantly smaller than 50 meV and so we show, in 
Fig. [2] the effect of reducing the electronic coupling. As V 
gets progressively smaller, the timescale of relaxation grows 
significantly, reaching approximately 100 ps for V = 1 meV. 
At the smallest values of V, panels (a) and (b), the relaxation 
rate can be seen to follow the expected k cc V 2 golden rule. 
Consistent with their perturbative origins, the NIBA dynam- 
ics become quantitatively exact for vanishing V, whereas Red- 
field theory's qualitative accuracy is maintained throughout all 
panels. Interestingly, for this value of the reorganization en- 
ergy, the secular and Markov approximations to the time-local 
Redfield equation yield impressive quantitative accuracy for 
all values of V. 

As another important numerical test, we now consider the 
effect of a third state on the dynamics of singlet fission, where 
an initial state couples to a second state which is in turn cou- 
pled to a third. This configuration is clearly akin to the me- 
diated mechanism, with the three states S\, CT, and TT. In- 
terestingly, the quantum dynamics of such mediating systems 



(a) Sequential (b) Superexchange 
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FIG. 3. Population dynamics of the three-state model described in the text. 
While both Redfield theory and NIBA provide a reliable description of dy- 
namics in the two-step sequential regime, NIBA is qualitatively unable to 
describe the superexchange regime. Redfield equation dynamics employ the 
secular and Markov approximations and exact results are calculated with the 
HEOM approach. 



has precedent in the donor-bridge-acceptor complexes of pho- 
tosynthetic charge transfer. Over 15 years ago, Makri and 
coworkers performed numerically exact path integral simula- 
tions of a three-state model very similar to the one considered 
here and detailed the dynamical features of two previously 
proposed, but physically distinct transport mechanismsi^i^ 
The first mechanism is a sequential one, whereby the inter- 
mediate state becomes transiently populated in a scheme well- 
described by two-step kinetics. This mechanism dominates in 
energetic regimes satisfying E\ > E2 > £3. A second mecha- 
nism, evincing the superexchange phenomenon, employs vir- 
tual states of the intermediate which is never directly popu- 
lated, yielding overall single-step kinetics for the 1 — > 3 trans- 
fer. Superexchange applies when the intermediate state energy 
is much higher than the other two, E% » E\ > £3. 

In light of the similarity with mediated singlet fission (both 
qualitatively and quantitatively, see below), we consider it 
an essential criterion that any adopted quantum dynamics 
methodology be able to capture this effect. To make connec- 
tion with the singlet fission problem, we will henceforth con- 
sider the definite state labeling referred to above, namely S 1 , 
CT, and TT. For comparison, we adopt the same electronic 
parameters used in Ref.|94|and the same system-bath coupling 
used above, H e i- P h = ~Li=s u CT,TT 10 *Lk,i c U<lu(i\- Again the 
baths have an Ohmic spectral density with Lorentzian cutoff 
parametrized by A — 25 meV, Ml =150 meV, and T = 300 K. 
We wish to strongly emphasize that although this Hamiltonian 
(electronic parameters to follow) was originally parametrized 
based on photosynthetic protein data, the magnitude of the 
parameters is almost identical with those expected of singlet 
fission. 

First, we consider the sequential regime, for which the 
electronic Hamiltonian has =0, Eqt = -50 meV, and 
Ejj = -250 meV, with V$ u ct = 3 meV and Vqtjt - 17 
meV. In Fig. [3ja), we plot the population dynamics of the 
three states and find that both Redfield theory and NIBA 
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agree quantitatively with each other and with numerically ex- 
act HEOM dynamics. Clearly, in this sequential mechanism, 
S i populates CT, which in turn populates TT. 

We now turn to the superexchange regime, for which 
£s, =0, Ect - 250 meV, and Ejj - -80 meV, with 
Vs u ct = Vctjt = 30 meV; note that the immense barrier, 
E C t - E Sl = 250 meV « 10k B T completely prohibits ther- 
mal activation. Turning to the results in Fig. [3] the Redfield 
and NIBA dynamics strongly differ, and only Redfield theory 
gives results in good quantitative agreement with the exact dy- 
namics. In the superexchange limit, the intermediate CT state 
is never directly populated, and the kinetics follows a sim- 
ple S i — > TT dynamical scheme. Thus, in spite of the rela- 
tively small electronic coupling values, Vs u ct = Vctjt - 30 
meV, superexchange must be understood as a higher-order ef- 
fect. The effective electronic coupling from Si to TT due 
to CJ-mixing is V e ff w Vs u ctVct,tt KEs t - Ett)- Perform- 
ing second-order perturbation theory with this effective elec- 
tronic coupling thus gives a rate which is overall fourth-order, 
explaining why superexchange eludes the usual second-order 
treatment, such as that employed in NIBA. Redfield theory, 
on the other hand, is completely nonperturbative in the elec- 
tronic couplings, and is thus entirely capable of capturing this 
highly relevant phenomenon. Although we will have more 
to say about it in our future work, this simple model clearly 
refutes arguments that high-lying CT intermediates preclude 
efficient mediated singlet fission. 



B. Results for spectroscopy 

Lastly, we apply the Redfield formalism to the calculation 
of linear absorption spectroscopy, the formalism of which is 
described in App. [D] In this situation, the non-Markovian 
time-local variant is to be preferred, as it exactly solves the 
so-called pure-dephasing problem appropriate for the single- 
molecule absorption spectrum. Employing the methodology 
described there, we have calculated the absorption of a pen- 
tacene molecule in solution, which compares very favorably 
with the experimental spectrum, see Fig. [4] In particular, the 
phonon sidebands are accurately reproduced, even though this 
is a purely non-Markovian, multi-phonon signature. As such, 
this feature cannot be described by the Markovian Redfield 
theory. Comparisons such as this one provide essential bench- 
marks for the accurate parametrization of both electron and 
phonon degrees of freedom, as well as the interaction between 
them. 

In the presence of intermolecular interactions, the absorp- 
tion lineshape will be changed from that of Fig. |4]due to the 
electronic coupling to other excited states. To discuss these 
effects, let us introduce the dipole operator which, to a rea- 
sonable approximation, is given by 

fi = J]\(S iW »s >(S \+H.c, (22) 

m 

where H.c. denotes the Hermitian conjugate. This approxi- 
mation follows from the observation that the transition dipole 
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FIG. 4. Calculated and experimental absorption spectrum of a single pen- 
tacene molecule at T = 100 K in solution. The calculation parameters are 
E(S i) = 2.3 eV, n = 170 meV, and A = 120 meV (S = 0.7), and a homo- 
geneous broadening of 30 meV has been applied. Experimental spectrum is 
from Refill 

matrix element for a single excitation from molecule m to 
molecule i, 

Hi„ = -e ^<0|r„|>r;;,> = -e {H m \r\L t ) , (23) 

n 

is significantly larger for m = i (Frenkel excitations) than for 
m + i (CT excitations), due to spatial locality. Thus only the 
intramolecular Frenkel excitation has a non-negligible transi- 
tion dipole moment, and charge-transfer and multiple-exciton 
states do not absorb light. However, the expression for the 
absorption signal (see Eq. ( ID 11 1) requires the adiabatic transi- 
tion dipole moments, which follow from the above equation 
and the diagonalizing transformation, \a) = 2/ C?|z), as 

M ff = <0|/i|ar> * Yj C (S0n^r (24) 

(Si)„ 

It is crucially important to recognize that the extent to which 
the adiabatic state a is composed of diabatic intramolecular 
excitations (C? s ■. ) determines the strength with which it ab- 
sorbs; this is the phenomenon of intensity borrowing. There- 
fore, signatures of "dark" diabatic states, such as charge- 
transfer or multi-exciton states, will appear in absorption spec- 
tra if these states are strongly coupled to the "bright" in- 
tramolecular Frenkel excitations. 

This phenomenon has been previously addressed in the 
MEG literature within the context of the coherent superposi- 
tion approximation, wherein the authors concluded that cou- 
pling to multi-exciton states does not affect the total absorp- 
tion coefficient, a. 15 While we agree that the integrated ab- 
sorption coefficient is indeed only determined by the bright 
singlet excitons, 

Jr»oo r*oo 
d<jL>I((jS) ~ y \/j. a £ \ 2 I d(i)8(u) - cj a o) 
o ^ Jo 

=zw 2= zzzfe.r c ^^i 2 (25) 

= 2 2 ((5i)m|(5i) " >|w ' |2=iV|j " 5 ' 12 ' 

(Si),„(Si)„ 
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we argue that the spectral structure (peak positions and in- 
tensities) is surely affected by the coupling to multi-exciton 
states. For example, consider an artificial system composed 
of the ground state, |0), a single bright state, \S\), and a 
dark (multi-exciton) state, \D), with equal energies E(S i) = 
E(D) = E and mutual coupling (Si\H e i\D) = V. The dipole 
operator is simply fl = /ilSiXOI + H.c. In the uncoupled 
limit V = 0, the electronic Hamiltonian is already diago- 
nal, and Eq. ( ID2I ) gives the absorption lineshape as I(to) = 
\fi\ 2 6(a> — E), and the absorption coefficient a = |ju| 2 . However, 
for V + 0, the adiabatic eigenstates are the symmetric and 
antisymmetric linear combinations, \a) = (2) _1 ^ 2 (|5i) + \D)), 
\p) = (2)- 1/2 (|5 1>-|£>», with energies E a = E-V,E P = E+V. 
In this case, the absorption lineshape shows two weaker peaks, 

1 9 

1(a)) = -M [<H W -E + V)+8(a>-E-V)] (26) 

but the same total absorption coefficient, a = \fi\ 2 . In light of 
the preceding analysis, we propose that evidence for coupling 
to dark, perhaps multi-exciton, states should be observable in 
the linear absorption spectrum, perhaps contrary to standard 
intuition. 

As a related point, we recall that there have been propos- 
als to experimentally seek out real-time quantum beating as 
evidence of coupling to multi-exciton states, 35 even if, in the 
limit of strong coupling, the frequency of beating becomes too 
high for experimental resolution in the time domain. How- 
ever, the origin of spectral peaks discussed above is identical 
to that of quantum beating, namely the oscillation of quantum 
coherences at frequencies given by energy differences. More 
importantly, in the frequency domain, the peak splitting is pro- 
portional to the strength of the coupling and thus more easily 
observed for strong coupling. In this sense, linear absorption 
and real-time quantum beating should be seen as complemen- 
tary tools for the investigation of coupling to multi-exciton 
states: weak coupling yields negligibly split peaks that may 
be difficult to resolve, but produces slow oscillations in real 
time that should be easy to observe. The situation is reversed 
for strong coupling, where spectral measurements should be 
preferred. As a proviso, the real-time observation of quantum 
beating may be an artifact of the diabatic basis. Specifically, 
the adiabatic populations may show pure exponential relax- 
ation, but because the transformation back to the diabatic ba- 
sis mixes populations with oscillatory coherences, the diabatic 
populations appear to exhibit quantum beating. Thus the real- 
time detection of such beating in part depends on the basis 
which is probed experimentally. 

The examples discussed in this section on spectroscopy and 
the preceding one on population dynamics illustrate the utility 
of a Redfield approach to the description of organic singlet fis- 
sion systems. In particular, most materials of current interest 
for singlet fission lie in a regime where the ratio Au/Mlu < 1, 
largely due to the dominant coupling to high-frequency aro- 
matic carbon bond stretching. For the same reason, these sys- 
tems generically have bath relaxation times that make the sim- 
plifying Markov approximation a sensible one. Lastly, Red- 
field theory and its variants are non-perturbative in the elec- 
tronic states. Thus, they do not alter the underlying descrip- 



tion of the frozen electronic structure theory and are capable 
of treating higher-order effects such as superexchange. This 
last point will be significant in our discussion of singlet fis- 
sion in pentacene dimers. 



V. CONCLUSIONS 

To summarize, we have presented a self-contained micro- 
scopic theory of multi-exciton formation in the context of sin- 
glet fission. Our formalism emphasizes the difference in elec- 
tronic bases, diabatic and adiabatic, as applied to both theoret- 
ical development and experimental interpretation. Building on 
this electronic foundation, we have applied techniques from 
the theory of open quantum systems to describe the finite- 
temperature quantum dynamics of charge and energy transfer 
processes taking place in singlet fission materials. Specifi- 
cally, such processes are facilitated by phonon absorption and 
emission, which can be given a microscopic origin in terms of 
certain vibrational motions of molecules. 

We furthermore discussed various approximate quantum 
master equations for the reduced density matrix and found 
that while NIBA-like theories which are perturbative in the 
electronic couplings yield accurate dynamics for two-state 
systems in regimes expected to hold in organic systems un- 
dergoing singlet fission, their perturbative nature is exposed 
by higher-order processes, such as the superexchange mecha- 
nism. On the contrary, Redfield-like theories, perturbative in 
the system-bath coupling, yield reasonably accurate results in 
essentially all regimes of relevance for singlet fission. We ad- 
ditionally elucidated the importance of a theory for both the 
populations and coherences of the reduced density matrix, al- 
lowing for investigation of dynamics in different electronic 
bases as well as prediction of linear spectroscopies such as 
absorption. While more accurate quantum dynamics scheme 
are certainly worthy of consideration, we emphasize the effi- 
ciency of Redfield theory, which allows for a rapid, but thor- 
ough, investigation of parameter space in both small and large 
systems, which is ideal for a computational screening of effi- 
cient singlet fission target materials. 

In the small model systems considered here, we have found 
that both direct and mediated mechanisms are plausible path- 
ways to efficient singlet fission. For reasonable electronic 
Hamiltonian parameters, the phonon degrees of freedom facil- 
itate fission on picosecond timescales. In particular, we drew a 
potentially useful comparison with charge transport in photo- 
synthetic donor-bridge-acceptor systems in the context of CT- 
mediated fission, wherein both sequential and superexchange 
mechanisms should be considered possible. We emphasize 
that our proposal for a unified framework for the microscopic 
treatment of singlet fission in organic systems is based on ac- 
curacy, efficiency, and physical transparency. In particular, we 
have generalized existing techniques and used physical argu- 
ments and numerical benchmarks to marry them for the pur- 
pose of a microscopic and accurate treatment of fission in sys- 
tems ranging from dimers to crystals. The companion paper to 
this one begins this program in pentacene dimers while future 
work will consider large aggregates and bulk crystals. 
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Appendix A: Energies and couplings in truncated CI basis 



Using the minimal basis presented in Sec. Ill Bl here we give 
formulas for the diagonal and off-diagonal matrix elements of 
the exact electronic Hamiltonian operator, 

Hel = hijcl^Cjjr + - Yj VijucJ^cl^ci^Ck^ (Al) 

ij,cr ijkl,{Tcr' 

where the sums are over all molecular orbitals of the isolated 
molecules (the indices i, j, k, I now contain both the molecule 
and its orbital) and 



v = j 



(r) 



4/r) = (i\h\j), (A2) 



V ija = J d 3 ri J rf 3 r 2 0*(ri)0}(r 2 )r^V*(ri)^(r 2 ) s (ij\kl) 

(A3) 

are one- and two-electron integrals over spatial orbitals. Here 
and henceforth we employ atomic units. 

The diagonal matrix elements of the Hamiltonian opera- 
tor in the diabatic basis ("energies") are given by expressions 
available from known CI theory, 

E[S ] = 2(i\h\0 + 2(ij\ij) - (ij\ji) (A4) 



and 



i,jeS 

E [(5i) m ] = E[S ] +E g + 2K HmLm - J HmLm , 
E [C m A„] = E[S ] + E g + 2K HmLn - J HmLn , 

E [T m T„] = E[S ] + 2E g + J HmHn + J LmLn 

- Jl m H m - J h m H u - Jl„H,„ - Jl„H„ 

-(l/2)(K HmHii +K LiifLn ) 

+ [(6V3+5)/4] (K LmHn + K LiiH J 

-[(6V3-5) /4] (K LmHm + K LiMi ) . 



(A5) 
(A6) 



(A7) 



We have introduced the gap E g = si m - Sh„ and notation for 
direct Coulomb integrals, Jn = and exchange integrals, 

While one could in principle calculate all possible cou- 
plings between the previously introduced diabatic states, those 
couplings involving three or four molecules will be negligibly 
small due to the weak overlap of the localized molecular or- 
bitals. Specifically, we propose to neglect all three-center and 
higher two-electron integrals. This semi-empirical approxi- 
mation should not drastically affect the results. 

Introducing the spatial orbital matrix elements of the Fock 
operator, F, 

(i\F\j) = (i\h\j) + 2{ik\jk) - (014*), (A8) 

keS 

the remaining off-diagonal matrix elements ("couplings") can 
be evaluated to give 

(C m A n \H e i\(S i) m ) = (L m \F\L n ) 

+ 2yH m L m \L n H m ) — (H m E m \H m E n ) 



(A9) 

(C m A n \H e i\{S 1),,) = -(H m \F\H n ) 

(A10) 

+ 2{H n L n \L n H m ) — (H n L n \H m L,^) 
{C m A„\H e i\C n A m } — 2(H ln L m \L n H„) — {H m L m \H n L n ) (All) 
{C m A„\H e ,\T m T n ) = ^{(L m \F\H n ) 

\H n L„) — (L m H m \H n H m )} 

(A12) 

{(S i) m \H e i\(S i)„) = 2(H m L n \L m H n ) - (H m L„\H„L m ) (A13) 

{(Si) m \H el \ > = Jhj2{(L m L n \H n L m ) (A14) 

- (H,„H„\L n H m )}. (A15) 

We have neglected the coupling to the ground state, e.g. 
(C,„A„\H e i\So). Note that such terms do not strictly vanish 
as they do in traditional CIS (Brillouin's theorem) because 
the reference state is not the Hartee-Fock solution of the full 
molecular cluster. Nonetheless, because the energy gap be- 
tween ground and excited states is large, the renormalization 
of excited states due to coupling to the ground state should be 
negligible. Furthermore, although these terms could in prin- 
ciple facilitate non-radiative decay to the ground state, we as- 
sume a bottleneck for phonon emission prevents such events, 
justifying our neglect of such couplings. 



Appendix B: Derivation of the system-bath Hamiltonian 

We begin by considering the nuclear dependence of the dia- 
batic electronic state energies and couplings introduced above, 
i.e. Ui(Q) = {i\U(Q)\i) and Vy(Q) = (i\U(Q)\j)- To simplify 
notation, we will employ the Roman characters i and j to de- 
note diabatic states and Greek characters (a,/3, y, ...) to denote 
the adiabatic eigenstates. 
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In the limit where the ground diabatic state is a local min- 
imum, we may perform a second-order Taylor expansion of 
the potential in terms of the mass-weighted coordinates Q„ = 



M,,(Q„ - Q„ ) and employ a transformation to the normal 
modes that diagonalize the Hessian q k — Yjn=\ 2jx=i m « x Q",x, 
yielding 



Uo({q„}) ^Eo + ^J] 



2 2 



(Bl) 



where Eq = Uq({0}) and {co^} are the eigenvalues of the 
Hessian. 53 By promoting the nuclear coordinates to operators 
and re-inserting the kinetic energy operator, we arrive at the 
ground state diagonal matrix element of an electron-phonon 
Hamiltonian in normal-mode coordinates, 



<0|£ f „r|0> = E + J] 



PI 1 2.2 

2 + 2^* 



(B2) 



These coupling constants can also be incorporated into an off- 
diagonal spectral density 



Jij(t*>) = 7: y, —S(a) - u> k ) 
2 k ^ 



(B8) 



with a corresponding "reorganization" energy Atj = 
7T J Q duJij{co)l CO. 

Combining all of the above, we thus arrive at the desired 
system-bath-type Hamiltonian, Eqs.fTTl 



Appendix C: Explicit expressions for the Redf ield tensor 
elements 



The Redfield tensor is given explicitly as 



SI, 82 



The matrix elements of the total Hamiltonian in the higher- 
lying diabatic states thus additionally acquire a linear term 
describing the shift of the excited state potential energy sur- 



face minimum located at {q, ''}, 



(i\6 tot \i) = Uittqf}) + X T + 1^ ^ k ~ q *) 



Pk 1 2^2 



(B3) 



where the vertical energy is Ej = f/,({0}) = Ui({q. }) + An, 
the Holstein-like coupling constants are given byS 5 -^ c k ,i = 
-cojqf and the reorganization energy 5 -^ of state i is defined 



k k k 

It will be convenient to now define the so-called spectral den- 
sity of the phonons, which completely characterizes the har- 
monic environment intrinsic in the normal-mode decomposi- 
tion. The spectral density of state i is defined by 



Jii(u) = —S(to - co k ) 
2 *-i lOk 



(B5) 



Rapysit) - r^ ar (f) + T my {t) 



(CI) 



where 



£ dse-^\ft^ p \ S )i%; p \G)) ph 



(C2) 
(C3) 



are integrals of thermal bath correlation functions. We have 
introduced the notation 



and (. . . } ph = Tr ph {. . . e' 6 *' 1 * 1 ' /Z ph }. The calculation of the 
thermal bath correlation functions required for the Redfield 
relaxation tensor is straightforward for the harmonic baths de- 
rived above. Assuming uncorrelated bath fluctuations, one 
finds 

^ c k jc k ,j{q k j(t)q k ,j(0)} ph = 6ijCu(t), (C5) 
k , k ' 



from which the reorganization energy can be rewritten as An - 

i poo ""' " 

7T J da)Jii(LO)/LO. 

The coordinate dependence of the off-diagonal matrix el- 
ements in these normal modes can then also be evaluated to 
first order, 



(i\6to,\j) = Vij + J] 



c k ,ijq k , 



with the Peierls-like coupling constants given by 55 - 98 
dV u ({q n }) 



Ckjj 



dqk 



(B6) 



(B7) 



~^^CkjjC k ^ mn {q k jj{t)q k ^ mn {0))ph = (Sjmdjn + 6j„6j„)Cij(t), 



k.k' 



(C6) 

and {q k jj(t)q k \ m (0)) p i, = 0. The functions Cy(f) are given by 
the usual weak-coupling expressions 

Cij(t) = — J duJij(co) jcoth j cos(tLtf) - i sin(wf) j • 

(C7) 

Observe that in the Markovian limit, Eqs. (!C2| i- (lC3t reduce to 
ordinary one-sided Fourier transforms. 
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Appendix D: Calculation of spectroscopic observables 

Spectroscopic observables, such as absorption and emission 
lineshapes, can be straightforwardly calculated from RDM 
dynamics. For example, the e-polarized absorption lineshape, 
I € (u>), is known to be given by the Fourier transform of the 
dipole-dipole correlation function, D e (f). 53 ' 99 Within certain 
approximations, 53 the dipole-dipole correlation function can 
be written as D € (t) = 2 ff \Ma,e\ 2 Pao{t) with p a o(0) = 1, where 
/i ff e is the e component of the transition dipole moment from 
the ground state to the electronic adiabatic state a, yielding 



„ Jo 



dte ,u " PM (t). 



(Dl) 



As a simple example, consider neglecting the electron-phonon 
coupling, such that p a o(t) = exp(-/w ff of); then the absorption 
spectrum takes the familiar form 



I e (co) = V \fx a , e \ 2 Re dte iM e~ k 



(D2) 



which clearly neglects any broadening or multiphonon effects. 
In the general case, i.e. by propagating the reduced den- 
sity matrix p(0) = |ar)(0| under the Redfield equation, such 
environmental effects may be included, as will be demon- 
strated numerically in Sec. [IV] Renger and Marcus" proposed 
the analytically useful procedure whereby one keeps the full 
time dependence of the diagonal coherence dephasing tensor, 
T^ aaa (f), and performs the Markov approximation for the off- 
diagonal tensors, r*^, all within the secular approximation, 
i.e. 



dp a o(t) 
dt 



-iu M p M {t) + R a oao(t)Pao(t), (D3) 



with 



RaOaoQ) * -T + mma (t) - £ Y + ma (t oo). (D4) 



Before concluding, we consider a specific example of a 
spectroscopic calculation within the time-dependent Redfield 
formalism, relevant to the absorption of a single molecule, the 
so-called pure dephasing problem . 100 ' 101 In this case, there is 
no electronic coupling, Vy = 0, and the Hamiltonian is com- 
pletely specified by the two level system, 



H tol = \S )E(S )(So\ + \Si) 



Ckqk 



(D5) 



VI 2 



2 2 



The time-local Redfield equation of motion for the coherence 
P5,s ( f )is simply 



dps, s {t) 
dt 



-i[E(SO-E{S y\p Sl So(f) 

+ RSiSoSrfoi^PSiSoit)' 



(D6) 





1.0 n 




- 

0.8- 








0.6 - 


>> 




= 


0.4- 




s 






0.2- 




o.o- 



(a) S = 0.5 



-r 



1.0 

0.8 
0.6 
0.4 
0.2 
0.0 



(b) S = 2.0 



1 2 3 4 5 

(ill - U1S 1 S + A Si) 



1 2 3 4 5 

(u-uiSiS a +A Sl )/n 



FIG. 5. Single molecule absorption spectra at T = for a single vibrational 
frequency, fl, with the Huang-Rhys parameter, 5 , as given. Spectra have been 
artificially broadened for clarity and normalized to the value of the 5 = 0.5 
zero-phonon (0-0) line. 



which can be straightforwardly solved. With the initial condi- 
tion p(0) = \S i)(Sq|, one finds 



PS,S (0 = ex P [-KWtSo - Asjt + g(t) - g(0)] , 
where the lineshape function g(t) is given by 
Js^oj) f , (/3tj\ 



(D7) 



g(t) 



1 f°° J s .(ai) ( lPto\ 
— — J da> j — jcoth y— J cos(wf) - i sin(wf) 



(D8) 

In fact, the pure dephasing problem can be straightfor- 
wardly solved exactly, 101 using the fact that the Hamiltonian 
is already diagonal in the electronic states and well-known 
thermal properties of harmonic oscillators. If one carries out 
this procedure, it is found to give exactly the same result as 
Eq. dD7t , a remarkable property of the time-local Redfield 
equation. Thus, the cumulant resummation inherent in the 
time-local formalism is exact for the pure dephasing problem. 
This result does not hold for the the time-nonlocal approach. 

To make our example more specific, consider a single vi- 
brational mode at frequency O, J(a>) = nS Q. 2 6(a> - Q), such 
that 

g(t) = S [coth(y6Q/2) cos(Qf) - i sin(fif)] • (D9) 

We have introduced the dimensionless Huang-Rhys factor, 
S - n 1 J Q da)J(a>)/to 2 making clear that g(0) is, at zero tem- 
perature, identical to the Huang-Rhys factor. Although an an- 
alytical evaluation of the required Fourier integral, Eq. (IDU . 
is still difficult, it can be straightforwardly calculated numer- 
ically. We show in Fig. [5] example absorption spectra calcu- 
lated for two different values of the Huang-Rhys parameter, S . 
Clearly the well-known vibronic progression is perfectly cap- 
tured, even in regimes of very strong system-bath coupling. 
Again we emphasize that this nonperturbative multi-phonon 
behavior is a purely non-Markovian effect which is only cap- 
tured exactly by the time-local form of the Redfield equation. 
The Markovian limit would yield only a single Lorentzian 
lineshape at a> — cos^o + dsCs^s), with broadening 
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